Temperature and Length Scale Dependence of Solvophobic Solvation in a Single-site 

Water-like Liquid 



John R. Dowdle 

The Dow Chemical Company, Freeport TX, 77541, USA 

Sergey V. Buldyrev 
Department of Physics, Yeshiva University, New York, NY 10033 USA 

H. Eugene Stanley 

Center for Polymer Studies and Department of Phystcs, Boston University, Boston, MA 02215 USA 

Pablo G. Debenedetti 

Department of Chemical and Biological Engineering, 
Princeton University, Princeton, New Jersey, 08544 USA 

Peter J. Rossky 

Institute for Computational Engineering and Sciences and Department of Chemical Engineering, 
The University of Texas at Austin, Austin, Texas 78712, USA 

The temperature and length scale dependence of solvation properties of spherical hard solvophobic 
solutes is investigated in the Jagla liquid, a simple liquid that consists of particles interacting via 
a spherically symmetric potential combining a hard core repulsion and a longer ranged soft core 
interaction, yet exhibits water-like anomalies. The results are compared with equivalent calculations 
for a model of a typical atomic liquid, the Lennard- Jones (LJ) potential, and with predictions for 
hydrophobic solvation in water using the cavity equation of state and the extended simple point 
charge (SPC/E) model. We find that the Jagla liquid captures the qualitative thermodynamic 
behavior of hydrophobic hydration as a function of temperature for both small and large length 
scale solutes. In particular, for both the Jagla liquid and water, we observe temperature-dependent 
enthalpy and entropy of solvation for all solute sizes as well as a negative solvation entropy for 
sufficiently small solutes at low temperature. This feature of water-like solvation is distinct from the 
strictly positive and temperature independent enthalpy and entropy of cavity solvation observed in 
the Lennard- Jones fluid. The results suggest that, compared to a simple liquid, it is the presence of 
a second thermally accessible repulsive energy scale, acting to increasingly favor larger separations 
for decreasing temperature, that is the essential characteristic of a liquid that favors low-density, 
open structures and models hydrophobic hydration, and that it is the presence of this second energy 
scale that leads to the similarity in the behavior of water and the Jagla liquid. In addition the Jagla 
liquid dewets surfaces of large radii of curvature less readily than the Lennard- Jones liquid, reflecting 
a greater flexibility or elasticity in the Jagla liquid structure than that of a typical liquid, a behavior 
also similar to that of water's hydrogen bonding network. The implications of the temperature and 
length scale dependence of solvation free energies in water-like liquids are explored with a simple 
model for the aggregation of solvophobic solutes. We show how aggregate stability depends upon 
the size of the aggregate and the size of its constituent solutes, and we relate this dependence to 
cold-induced destabilization phenomena such as the cold-induced denaturation of proteins. 



I. INTRODUCTION 



Among the many anomalous properties of liquid water 
is the solvation behavior of small apolar solutes, which 
is characterized at ambient conditions by an unfavorable 
entropy of transfer from vapor phase to water and an 
atypical decrease in solubility with increasing tempera- 
ture. This behavior contrasts with typical solvents, which 
more readily accommodate apolar compounds as thermal 
fluctuations increase. The enthalpy of transfer for non- 
polar solutes to low-temperature water is actually neg- 
ative and favorable, but the solubility is dominated by 
the entropic penalty. These characteristics change as a 
function of temperature and solute size. At sufficiently 



high temperatures the enthalpy is large and unfavorable 
and is only partially compensated for by favorable trans- 
fer entropies. Similarly, for sufficiently large solutes, the 
poor solubility is dominated by the unfavorable enthalpy 
associated with the formation of an interface, which over- 
comes the favorable entropy gain [T]. 

Recent theoretical work in the field of hydropho- 
bic solvation [2HZ] has refocused attention on the size- 
dependence of solvation free energy for small and large 
solutes, which is generally accepted to play a poten- 
tially important role in the formation and stabilization 
of many biological structures including proteins and cell 
membranes. Specifically, it was demonstrated that the 
solvation free energies of simple hard sphere solutes in 
water at ambient conditions undergo a crossover in size 
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dependence at about 1 nm [TJ . For solutes of size smaller 
than 1 nm, the solvation free energy scales with the vol- 
ume of the solute, while for larger solutes it scales with 
the surface area. This crossover behavior is general to all 
liquids far from the critical point and near liquid-vapor 
coexistence, but the length scale of the crossover in wa- 
ter is greater than that of simple liquids, such as a simple 
Lennard- Jones (LJ) liquid This longer crossover dis- 
tance is attributed to water's propensity to create avail- 
able space throughout its hydrogen bonding network. 

Traditional explanations of hydrophobic behavior, and 
water-like anomalies in general, place emphasis on the 
orientational interactions of water molecules (hydrogen 
bonding) and the accompanying tendency for tetrahe- 
dral structure. However, it has been demonstrated [8HTU] 
that water-like thermodynamic and structural anomalies 
can also be manifested by a recently introduced family of 
spherically symmetric potentials which possess two char- 
acteristic length scales (the Jagla model (HKH]), a hard 
core and a longer ranged soft core repulsion. Further, 
the Jagla model has also been shown to exhibit water- 
like solvation thermodynamics |13j . In particular, the 
solubility of simple hard sphere solutes in the Jagla liq- 
uid is a non-monotonic function of the temperature, and 
furthermore, a polymer composed of such hard spheres 
exhibits a solvent-induced collapsed state with a stability 
diagram in the pressure-temperature plane reminiscent 
of that of a typical globular protein in water [T3lU5| . 
These results confirm that orientational interactions are 
not necessary to produce these features of water-like sol- 
vation behavior [TMT%] and suggest that the presence of 
two competing length scales is a fundamental physical 
feature of hydrophobic hydration. 

Questions still remain, however, about the similarities 
between solvation in the Jagla liquid and water. In par- 
ticular, what are the energetic and entropic contributions 
to the solvation free energy in the Jagla liquid and are 
they similar to those of water? Over what length scales 
do the analogies in solvation behavior between the two 
liquids extend? Is the length scale crossover behavior in 
the Jagla liquid similar to that of other simple liquids, or 
does it also mimic that of water? In the present study, 
we address all of these questions using extensive Monte 
Carlo (MC) simulations of the Jagla liquid. In addition, 
we compare results for water and the Jagla liquid to re- 
sults for the LJ liquid wherever possible. In doing so we 
clarify what is indeed unique to water-like solvation and 
what is common in typical liquids. 

This paper is organized as follows. In Section [TTJ we 
describe the theoretical and computational methods used 
to calculate the thermodynamic quantities of interest. In 
Section [HI] we describe the interparticle potentials used 
and the details of the simulation protocols. The results of 
the calculations are presented and discussed in Sec. |IV[ 
and conclusions and future directions are given in Sec. 
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II. THEORETICAL & COMPUTATIONAL 
METHODS 

All solvation properties of a solute may be obtained 
once the excess chemical potential is known. Thus, our 
calculations focus on the evaluation of the excess chem- 
ical potential of a cavity solute, p x c , which is formally 
given by 



(1) 



where T is the temperature, k B is Boltzmann's con- 
stant, and Po(R) is the probability of finding a cavity 
of size R or larger around a randomly located point in 
solution. For sufficiently small cavities, po(R) may be 
evaluated directly via the test particle insertion method 
[T§1 |2"U] . In dense liquids, however, the probability of 
observing density fluctuations extreme enough to accom- 
modate cavities much larger than the solvent particles is 
exceedingly small, and test particle insertion is known to 
fail in this case [2Tj . 

There are several methods available for the evaluation 
of chemical potentials for large cavities (see e.g., 22 ), 
but for the Jagla and LJ fluids in this study we choose to 
use the revised scaled particle theory (RSPT) of Ash- 
baugh and Pratt 121] ■ Here we give only a brief 
overview of RSPT which closely follows that given in Ref. 
[25] . For more detailed descriptions the reader is referred 
toRefs. [23j[24]. 

RSPT improves upon classical scaled particle theory 
(SPT) [H[27] by including multi-body correlations. The 
essential idea behind both RSPT and SPT is that the 
excess chemical potential must be equal to the work re- 
quired to inflate a cavity against the solvent from size 
zero to R. This work must oppose the pressure due to the 
solvent molecules at the cavity boundary, and thus scaled 
particle theories require knowledge of the contact corre- 
lation function, G(R), defined to be the average density 
of solvent molecules, relative to the bulk, at the cavity- 
solvent interface. With G(R) known, the excess chemical 
potential is calculated as 



(2) 



H X C {R) = / k B T P G(r)4Trr z dr, 



where p is the bulk solvent number density. For R 
much greater than the solvent size, the contact correla- 
tion function may be expanded in curvature, R , with 
phenomenological coefficients 
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(3) 



Here, P is the bulk pressure, is the surface tension 
of a flat solvent-cavity interface, and S is the first-order 
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curvature correction to the surface tension [53]. An ex- 
pression for the excess chemical potential of large cav- 
ity solutes is then be obtained by expanding Eq. Q to 
fourth order and integrating to get 



Mc(^)llargo = g + 47Ti? 2 7oo - 167T 7oo <5i? 

+ AirkgTpK , (4) 

H 

where A is the fourth-order curvature correction coeffi- 
cient and k is an integration constant. Third order coeffi- 



cients are typically set to zero so as to avoid logarithmic 
contributions to p x c [2H1 [29], a convention we follow in 
this work. The results for the test particle insertion cal- 
culations for small cavities, /z^(i?)| S i m , are interpolated 
with the large cavity solute expression in Eq. Q by 



(4(R) = p X (R)\ sin J(R) + Mc(fl)|large(l - /(#))■ (5) 

The function f(R) used here is a cubic function de- 
signed to smoothly switch between small (i? s im) and large 
(-Rlarge) cavity sizes, 
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In order to obtain parameters appearing in the expan- 
sion for the contact correlation function, we use Eq. (§ 
and differentiate Eq. ^ with respect to R to obtain the 
contact correlation function as 
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and fit this function to the contact values calculated 
from the MC simulations, as demonstrated in Fig. [T] 
The pressure is set equal to the simulation pressure, and 
the parameters 700, 5, k, and A are fit to the simulation 
results. 

The contact density calculations for the Jagla and LJ 
liquids demand significant amounts of computer time to 
obtain good statistics, and performing similar calcula- 
tions for typical multi-site water models that have elec- 
trostatic interactions is not desirable. For our purposes 
of comparison here, we may, however, estimate the excess 
chemical potential of large cavities in water over a broad 
range of thermodynamic states by using the recently de- 
veloped cavity equation of state (C-EoS) [3D]. The C- 
EoS is an analytical equation of state parameterized to 
fit experimental and simulation results for water, and 
it has been shown to accurately reproduce hydrophobic 
solvation thermodynamics of simple hydrophobes when 
combined with a first-order perturbation theory. The 
functional form of the C-EoS is given by 




Figure 1. Demonstration of a fit of Eq. |7| for the cavity 
contact correlation function to calculated contact values for 
several cavity sizes in the Jagla liquid at T = 0.6 [£2 /As]. The 
contact correlation function, G(R) (dashed line), is fit to the 
maxima (open circles) in the cavity-solvent pair correlation 
functions, gns-jGir) (solid lines). The cavity radii are, in 
units of o-jg, 0.32, 0.71, 1.09, 1.47, 1.86, 2.24, and 2.63. 



Pul = a + bp + cln/3, 



(8) 



where p, x is the cavity chemical potential and the co- 
efficients, a, b, and c are assumed to be temperature in- 
dependent. Thus, the C-EoS assumes that the enthalpy 
of cavity formation depends linearly upon temperature 
and that the associated heat capacity is temperature in- 
dependent. The dependence of \x x on the cavity size, R, 
is obtained by expanding in powers of 1/R and requiring 
that (3p x approach 7 ;„ao in the large cavity limit, where 
jl v is the experimental liquid-vapor surface tension and 
ao = 4irR 2 is the cavity surface area, 
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The remaining coefficients Ai, Bi, and Ci are obtained 
from fits to simulation data. 



III. SIMULATION DETAILS 

MC simulations of cavity solvation in the Jagla and LJ 
fluids were performed along the liquid vapor coexistence 
curves of each liquid for states ranging from the triple 
point to slightly below the critical point. The Jagla po- 
tential is given by 



where 



oo, r < r , 
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m 2 r + 6 2 , T\ < r < r 2 , 

0, r > r 2 , 
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This potential, shown in Fig. [2j demonstrates a wide 
range of behavior for varying parameters, including lim- 
iting cases of hard sphere, triangle well, and ramp po- 
tentials. Here we choose r\ = 1.72r , r 2 = 3.0r , and 
El = 3.5e 2 , as this particular parameterization manifests 
a cascade of water-like anomalies [9] H31 [3lJ [32] ■ 

For the LJ fluid we use the cut-shifted LJ interaction 
given by 



u LJ (r) - u LJ (r c ), r<r c 
0, r > r c 



(15) 



where u LJ (r) = 4e LJ (ofj/r 12 - crlj/r 6 ) is the full LJ 
interaction, elj and <jlj are the well depth and solvent 
diameter, respectively, and the cutoff distance, r c , used 
is chosen as 2.5<7lj. 

Several different sets of Monte Carlo simulations were 
performed on the Jagla liquid. In the first, saturation 
properties of the Jagla fluid were estimated from canon- 
ical ensemble MC simulations of a liquid slab in equi- 
librium with its vapor for selected temperatures ranging 
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Figure 2. The Jagla two-ramp potential. The parameters 
used in the present studies are the same as in [T3], viz.: r\ = 
1.72ro, r2 — 3.0ro, and ei = 3.5e2- The relative values of 
the hard core (ro) and the soft core (ri) positions roughly 
correspond to the same ratio between the positions of the 
first and second solvation shells of liquid water. The effective 
size of the Jagla particle, ujg, is estimated from plots of the 
radial distribution to be the minimum separation at which 
ujg(t-) = (see Fig. [3). 



from near the triple point to just below the critical point. 
From these slab simulations we estimate saturated liq- 
uid and vapor densities, the saturation pressure, and the 
liquid-vapor surface tension along the binodal line. The 
surface tension, 7^ , is calculated from the profiles of the 
pressure tensor using the mechanical definition [331 Hi] • 
The results for the saturation properties are shown in 
Table H 

In the second set of simulations, isothermal-isobaric 
MC simulations of the Jagla fluid were performed for 
both the liquid and vapor phases at each of the satura- 
tion states listed in Table SI in the supplementary mate- 
rial. Test particle insertion calculations were performed 
on the resulting liquid phase trajectories for cavities up 
to 2<t jq in diameter to obtain fJ.^(R)\ s i m . Similarly, in- 
sertion probabilities and excess chemical potentials for 
cavities up to 6a jq in diameter were obtained from test 
particle insertion analysis of the vapor phase trajecto- 
ries. Knowledge of the vapor phase chemical potentials 
allows evaluation of the surface tension at the vapor wall 
interface [25] . 

Finally, isothermal-isobaric MC simulations of a single 
cavity in the Jagla liquid were performed for various cav- 
ity sizes at each of the saturation states listed in Table SI 
in the supplementary material. Cavity diameters up to 
6(7/07 were considered, and the contact correlation func- 
tion was evaluated for each cavity at each state point. 
The contact correlation function is determined by ex- 
trapolating the cavity-solvent pair correlation function 
to contact. 

The parameters in Eq. Q may be fit to the MC re- 
sults for G(R), and the cavity excess chemical potential 
may then be computed from Eq. The details of the 
MC simulations used to calculate the insertion probabil- 
ities and contact correlation functions in the Jagla fluid 
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are provided in Tables SI and S2 in the supplementary 
material. All data for the LJ liquid are those obtained in 
the studies reported in Ref. [25] ■ The saturation states 
for the LJ liquid are also listed in Table S4 in the sup- 
plementary material for the present study. 

Molecular dynamics simulations of the SPC/E water 
model |35] were performed along the liquid vapor coex- 
istence curve for each of the states listed in Table S5 
in the supplementary material. A system consisting of 
512 SPC/E water molecules was simulated in a cubic 
box with periodic boundary conditions in the canonical 
ensemble for 20 ns using the GROMACS molecular dy- 
namics engine (36] [37] . The time step was chosen as 2 
fs, and bonds were constrained with the SETTLE algo- 
rithm }38| . The velocity rescaling thermostat was used to 
control temperature with a time constant of 0.1 ps [39] - 
Particle mesh Ewald summation was used to treat long 
range electrostatic interactions [ID] with a real space cut- 
off of 1.2 nm and a mesh spacing of 0.18 nm. The Ewald 
tolerance was set to 10~ 5 , and fourth order interpolation 
was used. 



IV. RESULTS AND DISCUSSION 

A. The Definition of Solvent Size from Pair 
Distribution Functions 

A comparison of the solvent-solvent pair correlation 
function, g(r), for the three liquids is shown in Fig. [3] 
The maximum in g(r) for the LJ liquid occurs at a pair 
separation slightly larger than a^j, and at a separation 
of cr^j the pair distribution function assumes a value of 
very nearly one for all states on the saturation curve. 
The nearest separation at which g(r) is unity is a com- 
monly used estimate for the size of a particle since the 
surrounding fluid is depleted from all shorter distances. 
We adopt this estimate here and use <7l,j as the size of 
the LJ particle. 

In the case of SPC/E water, the pair distribution func- 
tion peaks at about 0.28 nm at ambient temperature and 
slightly larger distances at higher temperatures. These 
distances are smaller than the LJ diameter for oxygen 
due to H-bonding. The nearest separation at which g(r) 
is unity is nearly constant at about 0.26 nm, which, to 
be consistent, is our choice for the size of the SPC/E 
molecule, aw at- 

The maximum peak in the Jagla liquid g(r) occurs at 
a distance significantly larger than the hard core diam- 
eter, tq. This reflects the preference of Jagla particles 
to maintain separation at the minimum in ujQ(r), n, 
unless stressed by temperature or pressure. This prefer- 
ence is diminished as temperature increases. However, 
the minimum separation at which the Jagla g(r) is unity 
is found to be insensitive to temperature [see Fig. [3] (c)] 
and closely corresponds to the minimum separation at 
which the pair potential is zero. This distance, ctjg, is 
a consistent estimate for the size of the Jagla particle; 
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Figure 3. Solvent-solvent pair distribution functions for states 
along the saturation curves of (a) the LJ liquid, (b) SPC/E 
water, and (c) the Jagla liquid. It is evident from the figure 
that the minimum separation at which g(r) has the value 
unity can be used as an estimate for the solvent size. For the 
SPC/E model this corresponds to aw at = 0.26 nm, for the 
LJ liquid it is <Tlj, and for the Jagla liquid it is ojq — 1.56ro 
(the minimum separation at which ujc(r) = 0). These sizes 
are taken to be independent of temperature for the states 
considered here, as justified by the data shown. 



<7jg = 1.56ro for the potential parameterization consid- 
ered here. 



B. Surface Tension and Vapor-Liquid Equilibria in 
the Jagla fluid 

In the first set of MC simulations, saturation prop- 
erties of the Jagla fluid were estimated from canonical 
ensemble MC simulations of a liquid slab in equilibrium 
with its vapor for selected temperatures ranging from 



G 



near the triple point to below the critical point. From 
these simulations we estimate liquid and vapor densities, 
the saturation pressure, and the liquid-vapor surface ten- 
sion along the binodal line. 

The results for the liquid-vapor slab simulations of the 
Jagla fluid are summarized in Table [I] The saturated liq- 
uid densities and the equilibrium vapor densities are in 
close agreement with those reported by Lomba et al. [H] . 
We expect that our estimates of the coexistence proper- 
ties of the Jagla fluid may be improved upon by taking 
finite size effects into account, as it is known, e.g., that 
large wavelength fluctuations may be suppressed by the 
system size [32] • Nevertheless, the solvation behavior we 
seek to characterize occurs for states at or near coexis- 
tence [3J, and we therefore expect the present estimates 
from the slab simulations to suffice for this study. 
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Table I. Canonical ensemble MC simulations of a liquid slab in 
equilibrium with its vapor were performed to obtain estimates 
of saturation properties. N Jagla particles were simulated at 
five different temperatures for 1.6 x 10 6 MC cycles, where one 
cycle corresponds to N MC moves. The liquid and vapor den- 
sities were estimated from ensemble averages of the densities 
in the centers of the liquid and vapor regions, respectively. 
Similarly, the saturation pressure was obtained by evaluating 
the pressure tensor in the center of the vapor region. The 
liquid-vapor surface tension is calculated using the virial re- 
lation [331 134| . Numbers in parentheses are estimates of the 
statistical error in the last digit of the reported value. 



C. Cavity Solvation Thermodynamics 

The parameters in Eq. ^ were fit to the MC results 
for G(R) in the Jagla liquid the using a least-squares re- 
gression. The choice of i? s i m and i?i a rgo used in the fit var- 
ied with the thermodyanmic state. Values of Ri m ranged 
from 0.5 to 0.6a jg and values of i?i ar go ranged from 0.75 
to 0.95a jg- In all cases, G(R) was well represented be- 
tween R lm and i?i ar gc by differentiation of /ij? (i?)| S i m . 
The results of the fit are presented in Table [TT] The sur- 
face tension of the flat interface, 700, is higher than the 
liquid-vapor surface tension measured in the slab simula- 
tions at all temperatures. It should be emphasized that 
7oo does not strictly correspond to the liquid-vapor sur- 
face tension, but rather to the total interfacial free energy 
between the solvent and the cavity which consists of con- 
tributions from two interfaces — a liquid vapor interface 
between the solvent and vapor film surrounding the cav- 
ity and the vapor-wall interface between the vapor film 



and the cavity surface. If the two interfaces are well sep- 
arated and not interacting with one another, then 7^ is 
equal to the sum of the liquid- vapor and vapor-wall sur- 
face tensions. Our simulations are sufficiently far from 
the critical point that the vapor-wall surface tensions are 
negligible for all states considered. Furthermore, the fit- 
ted values of 700 were insensitive to varying the maxi- 
mum cavity diameter used in the fits between 4erjG and 
6a jg, suggesting the finite cavity sizes considered here 
are not to blame. Therefore the difference between 700 
and ji v is likely due to other factors such as the finite- 
size limitations of our estimates of or the physical 
impact of quenched fluctuations at the solvent-wall in- 
terface [25]. The first order curvature correction, S, is 
negative and decreases with increasing temperature, also 
consistent with the results for the LJ liquid. It should 
also be mentioned that here S need not correspond to 
the Tolman length [43], but rather is treated as a fitting 
parameter. The parameters k and A are negative for all 
states and diminish in magnitude as the critical point is 
approached. 
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Table II. Parameters from the least-squares fit of Eq. |7| 
to the contact densities obtained from the simulations in Ta- 
ble S2 in the supplementary material. The simulation data 
was split into several blocks, and the numbers in parentheses 
represent an error in the last digit in the fitted parameter cor- 
responding to one standard deviation of the block averages. 

The results of the MC calculations for the cavity con- 
tact correlation functions are shown in Fig. [4] along with 
the fits to G(R). In both fluids, as the solute size grows 
from zero, the solvent packs increasingly tightly until the 
contact density peaks at a value of R on the order of 
the solvent size. At this point, the solvent begins to 
pull away from the solute, and for sufficiently large so- 
lutes, G(R) will be less than one. The contact corre- 
lation function is a decreasing function of temperature 
for all solute sizes studied here, but for sufficiently large 
solute sizes G{R) will increase with temperature since 
lmiR_ s . 00 G(R) = fiP/ p, which increases with tempera- 
ture along the saturation curve. 

The cavity sizes where G(R) decreases below one, i.e. 
where the cavity is "dewet", are larger relative to the 
solvent size in the Jagla liquid, meaning that the Jagla 
liquid resists dewetting of hard surfaces more than the 
LJ liquid. Lastly, for a fixed cavity size in the LJ liquid 
the spacing in G(R) values between temperatures ap- 
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(a) Lennard- Jones 




Figure 4. Cavity contact correlation functions as a function 
of cavity size (measured in units of solvent diameters) for 
states along the saturation curves of the (a) LJ and (b) Jagla 
liquids ranging from near the triple point (blue) to just below 
the critical point (red). The temperatures for the LJ liquid 
range from ksT/sLj = 0.65 (blue) to 1.00 (red) in increments 
of 0.05, while those for the Jagla liquid range from ksT '/e2 
= 0.4 (blue) to 1.2 (red) in increments of 0.1. Points are 
obtained from MC simulation data and lines are fits of Eq. 
Q to the simulation data. Statistical errors are smaller than 
symbol size. All LJ data are obtained from Ref. |25| . 



pears roughly constant, suggesting a linear dependence 
upon temperature. This is not the case in the Jagla liq- 
uid, however, as the temperature dependence clearly de- 
creases with increasing temperature. 

With the fitted parameters for G(R), the excess chem- 
ical potentials for the Jagla and LJ liquids may be ob- 
tained from Eq. In the case of water we use Eq. 
([9]). The results of the chemical potential calculations 
are shown in Fig. [5] The excess chemical potential is a 
positive, monotonically increasing function of cavity size 
at all temperatures in all three liquids. 

In the LJ liquid, the chemical potential is a decreasing 
function of temperature for all cavity sizes greater than 
<tlj/2. Furthermore, the spacing between temperatures 
for any fixed cavity size appears roughly constant in the 
LJ liquid, which, as pointed out by Ashbaugh [25], sug- 
gests that along the saturation curve the excess chemical 
potential may be modeled as 



fi(R) = h* D (R)\ < ,-T&R)\„, 



(16) 
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where h x c {R)\ a and s x c {R)\ a are the temperature inde- 
pendent enthalpy and entropy of solvation. The enthalpy 



R 



Figure 5. Excess chemical potential per surface area versus 
cavity size (measured in units of solvent diameters) for states 
along the saturation curves of (a) the LJ liquid, (b) water, and 
(c) the Jagla liquid. The thermodynamic states for the LJ and 
Jagla liquids are the same as those presented in Fig. [4] Points 
in the Jagla and LJ plots are obtained from simulation data 
and scaled particle theory. Lines in the LJ plot are fits using 
Eq. ( 16 1, while lines in the Jagla plot are fits of the simulation 
data to the C-EoS [Eq. (JoJ) ] . Lines in (b) are predictions from 
the water C-EoS [3U]. The temperatures used for the water 
C-EoS plot are T [K] = 273, 304, 335, 366, 398, 429, 460, 
491, and 522. 



is positive and increases with cavity size, indicating the 
loss of favorable solvent-solvent interactions near the cav- 
ity solute. Except for cavities smaller than <Jlj/2, the 
entropy is also a positive, increasing function of cavity 
size, indicating that solvent molecules near the cavity 
experience a net gain in configurational space. The ex- 
cellent fit of Eq. (16) to the simulation data [Fig. |5ja)], 
indicates that that the enthalpy of solvation is approxi- 
mately temperature- independent, and therefore the heat 
capacity of cavity solvation in the LJ liquid is approxi- 
mately zero. In the Jagla liquid, in contrast, the chemi- 
cal potential is an increasing function of temperature for 
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small, solvent-sized cavities and a decreasing function of 
temperature for large cavities. The temperature deriva- 
tive of the excess chemical potential for a fixed cavity 
size is not constant [Fig. [jjjjb)], but is evidently nonlin- 
ear. The qualitative behavior of the chemical potential 
of cavity solvation in the Jagla liquid is remarkably simi- 
lar to that predicted for liquid water by the C-EoS. This 
suggests that the Jagla liquid data may be fit to the C- 
EoS as well. Using the surface tension data (Table [I]) and 
a least-squares fit of the excess chemical potentials cal- 
culated from the G(R) data, we obtained a set of C-EoS 
parameters for the Jagla liquid (see Table S5). The fit 
is, in fact, excellent for all cavity sizes and temperatures 
considered, with slight deviations occurring only for the 
largest cavities at the highest temperature. The C-EoS 
fit to the simulation data permits exploration of the ther- 
modynamic contributions to [x x c in the Jagla liquid using 
analytical derivatives of Eq. {(oj) . 

The enthalpic and entropic contributions to the excess 
chemical potential for the Jagla liquid and water may be 
obtained from analytical temperature derivatives of the 
C-EoS Q The enthalpy and entropy of cavity solvation 
are compared in Fig. [6j The most obvious distinction be- 
tween the three liquids is that the LJ liquid has temper- 
ature independent enthalpic and entropic contributions 
to the solvation free energy, while the contributions for 
the Jagla liquid and water both show a strong tempera- 
ture dependence. For all three fluids, the enthalpy is a 
positive, monotonically increasing function of the cavity 
radius. The unfavorable enthalpy results from the disrup- 
tion of the liquid structure in the vicinity of the solute 
and the concomitant formation of an interface which on 
average has fewer favorable solvent-solvent interactions 
than an equivalent volume in the bulk. 

For any fixed cavity size in the size ranges considered 
in this study, the enthalpy is an increasing function of 
temperature in the Jagla liquid and in water. A pos- 
sible interpretation for this result in water is given by 
the Muller model [351 ES] > which uses a simple two-state 
hydrogen bond (H-bond) model parameterized by empir- 
ical solvation data to argue that the fraction of broken 



1 The temperature derivatives are taken along the saturation 
curve, cr, and they may be related to their constant pressure 
counterparts through the state variable relation 1441 



9}4 
dT 



d}4 
dT 



Noting that (d^/dP) T = 
molar volume, we may write 
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dP 



dp 



where is the excess partial 



BP 
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Similarly, 
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c \dT 



The fundamental differences between the liquids considered here 
are seen in both the saturation and constant pressure quantities. 



H-bonds in the solvation shell of apolar solutes is always 
at least somewhat greater than that in the bulk, and fur- 
thermore, that this disparity increases with temperature. 
Thus, for a fixed cavity size an increase in temperature 
decreases the number of H-bonds in the solvation shell 
relative to the bulk, which leads to a greater enthalpy. 

The entropy of cavity formation in both the Jagla liq- 
uid and water increases with increasing temperature for 
any fixed cavity size. It is possible that this behavior in 
water may be also be connected to the breaking of solva- 
tion shell H-bonds. If an increase in temperature causes a 
decrease in the number of solvation shell H-bonded pairs 
relative to bulk, then overall the gain in configurational 
freedom will be larger at the higher temperature. How- 
ever, this does not yet explain the Jagla model behavior. 

It is remarkable that the Jagla liquid, which contains 
no orientational dependence in its interaction poten- 
tial and therefore no H-bonding, reproduces the qualita- 
tive behavior of hydrophobic hydration thermodynamics. 
The underlying physical origins for this behavior in the 
Jagla liquid may be analogous to those of water, however. 
It has been shown in computer simulations of SPC/E wa- 
ter that the energetics of H-bonding are strongly corre- 
lated with local crowding effects. In particular, H-bonded 
pairs with a small number of neighbors will on average 
have a stronger H-bond than bonded pairs with a greater 
number of neighbors [47] . Furthermore, the fraction of 
H-bonded pairs in interfacial regions of apolar moieties is 
lower than in the bulk liquid, and the bonded pairs that 
do exist in these regions tend to have fewer neighbors 
and stronger bonds than the average H-bonded pair in 
the bulk. The interpretation is that density fluctuations 
that create cavities select against weak H-bonds, leaving 
only the stronger bonds to survive. Thus, the interfa- 
cial region experiences less H-bonding on the whole than 
equivalent volumes in the bulk, but maintains on average 
stronger hydrogen bonds. 

A plausible analogy in the Jagla liquid to H-bonding in 
water is the interaction of particle pairs at the potential 
minimum distance, r±. As temperature is lowered, the 
liquid prefers to adopt configurations that maximize the 
number of particle pairs near a separation of rx, which 
in the limit of the crystal is an hep lattice [32] • This is 
analogous to water maximizing the number of H-bonded 
pairs at low temperatures by adopting a tetrahedral net- 
work structure, and thus the Jagla pair interactions near 
r\ become analogous to water's H-bond. Under this view, 
density fluctuations in the Jagla liquid disrupt weakly in- 
teracting Jagla particles and leave a solvation shell that 
consists of fewer pair interactions near r\ . The fraction 
of "broken" interactions at r\ in the solvation shell would 
increase faster with temperature than the same quantity 
in the bulk. Future work entailing a detailed analysis of 
solvation shell structure will be needed to demonstrate if 
this hypothesis is correct. 

The LJ liquid demonstrates enthalpic and entropic be- 
haviors in sharp contrast to those of water and the Jagla 
liquid. The entropy is strictly positive for all cavities 
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Figure 6. (sl-s3) Entropy and (hl-h3) enthalpy of cavity solvation for the LJ liquid, water, and the Jagla liquid as a function 
of cavity size (measured in units of solvent diameters) . The temperatures for the Jagla liquid are the same as those listed in 
Fig. [4] while the temperatures for the water C-EoS are the same as those listed in Fig. [5] For water and the Jagla liquid, 
entropies are calculated from temperature derivatives of the cavity equation of state (lines), s^L = —(dfj^/dT)^, while for the 
LJ liquid, the entropy is given by the assumed temperature-independent form of fx® in Eq. ( ]16[ l. The enthalpy is calculated 
from hc\ a = + Ts^.\ a . Points in (s3) and (h3) are numerical derivatives of cubic spline fits to the excess chemical potentials 
in Fig. [5] 



of size R > 0.5ctlj in the LJ liquid, and the heat ca- 
pacity increment is negligible. This latter phenomenon 
is consistent with the argument for the temperature de- 
pendence of the relative fraction of broken H-bonds in 
solvation shell water compared to bulk water — i.e., the 
absence of a second energy scale in the LJ liquid pre- 
cludes a temperature-dependent enthalpy of cavity for- 
mation analogous to that of water. This implies that the 
fundamental commonality between water and Jagla flu- 
ids is the presence of two energy scales, each coupled to a 
different length scale, so that low density, open structures 
are increasingly favored for decreasing temperature, the 
feature absent in simple liquids. In the Jagla model, the 
second energy and length scale is set by the ratios de- 
scribing the soft ramp, £1/^2 and ri/rQ , while in water, 
these are determined by characteristics of the H-bonded 
and non- H-bonded states. 



D. The Length Scale Crossover 

As seen in Fig. [5j the chemical potential decreases with 
temperature along the coexistence curve for all cavity 
sizes considered in the LJ liquid. However, in water and 
the Jagla liquid, the chemical potential increases with 
increasing temperature for solvent-sized cavities and de- 
creases with temperature for larger cavities. Qualita- 
tively, the temperature dependence of the solvation free 
energy is identical in the Jagla liquid and water. 

An important consequence of the similarities between 
the temperature-dependence of the solvation free energies 
in the Jagla liquid and water is that the water-like char- 
acteristic of negative solvation entropy for small cavities 
is observed in the Jagla liquid (Fig. [6j). As the cavity 
size increases from R = 0.5a.jc, the curves along each 
saturation state first decrease, then pass through a mini- 
mum before increasing monotonically for larger cavities. 
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For cavities large enough that s^\ a > 0, the solvation 
shell is more disordered, and for sufficiently large cav- 
ities a dewetting transition will occur. This "entropic 
crossover" from negative to positive solvation entropy 
may therefore be viewed as a measure of the length scale 
at which interface formation begins to dominate the sol- 
vation free energy. In this view, the crossover for the 
LJ liquid occurs at cavity sizes less than a^j in diam- 
eter for all saturation states, which is smaller than the 
smallest cavities explicitly studied here. In water and 
the Jagla liquid however, the entropic crossover distance 
grows many times larger than the solvent diameter as 
temperature is decreased, as shown in Fig. [7| Although 
the entropic crossover is similar in the Jagla liquid and 
water, the crossover in water occurs at larger sizes rela- 
tive to the solvent diameter. 
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Figure 7. Entropic sign crossover lengths for cavity solutes in 
the Jagla liquid and water as predicted by the cavity equation 
of state. The values are plotted as a function of temperature 
reduced by T cr it, the liquid- vapor critical point. Points indi- 
cate the cavity radius, in units of solvent diameters, at which 
the solvation entropy changes sign from negative to positive. 
Entropic crossovers for cavities in the LJ liquid also occur, but 
at cavity radii less than 0.5 solvent diameters for all states on 
the saturation curve (not shown). 



AG = [IB, - nfi r , 



(17) 



where [Ir is the aggregate's chemical potential and \x r 
is the chemical potential of a single constituent solvo- 
phobe at infinite dilution. The relationship between the 
number of hard spheres comprising the aggregate and its 
radius, R, is n — Attt/R 3 /3v. Combining the expressions 
for n and AG and dividing by the aggregate surface area, 
we have 



AG(R) /AttR 2 = n R {R) /4irR 2 - ^R/Zv. (18) 
For increasing R, the first term on the RHS of Eq. 



(18) becomes approximately constant and equal to the 
interfacial free energy per unit area [3] . The second term 
is a linear function of the aggregate radius. The radius at 
which the RHS vanishes is the aggregation radius, R a — 
aggregates of size larger than R a are thermodynamically 
stable within this model free energy. These concepts are 
shown pictorially in Fig. [8] 
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E. The Thermodynamic Stability of Solvophobic 
Aggregates 

To explore the implications of the interplay between 
temperature and length scale dependence of solvation 
free energies, we examine a simple picture of solvophobic 
aggregation that combines ideas from Chandler [T] and 
Rajamani et al. [Hj. Consider a solvophobic aggregate 
composed of n identical hard sphere particles with cavity 
radius r such that the total volume of the aggregate is 
V = nv/rj, where v is the volume of a single constituent 
hard sphere particle and rj is the packing fraction of the 
spheres. If the aggregate is treated as a large spherical 
volume of radius i?, then the aggregation Gibbs energy 
may be modeled as 



Figure 8. Solvation free energy scaled by the surface area 
versus aggregate radius. The solid line correspond to the sol- 
vation free energy per unit surface area of a cavity of size 7?, 
which is used to model an aggregate of n smaller hard spheres 
of size r (see text). The dashed line represent the solvation 
free energy per unit surface area for n constituent spheres 
fully dispersed in solution. Only aggregates larger than the 
aggregation radius, R a , are thermodynamically stable. 

We now consider the process of cooling the aggregate 
from a warm temperature, Th, to a lower temperature, 
Tl, and in particular, the effect that this process has on 
the thermodynamic stability of the aggregate. A quali- 
tative picture of the dependence of the aggregation ra- 
dius, R a , on temperature for a water- like and a refer- 
ence LJ-like fluid is shown in Fig. [9] The differences in 
crossover behavior arise due to the fact that for small 
solutes in water-like solvents, increasing the temperature 
decreases the solubility. This has two effects: the first is 
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that the crossover length scale is more sensitive to tem- 
perature, and the second is that the slope of the dis- 
persed solvophobes line for high temperature is greater 
than the corresponding line at low temperature. These 
effects combine to produce a range of aggregate sizes that 
are thermodynamically stable at Tjj but become unsta- 
ble upon cooling a to Tj,. It is interesting that such a 
region also appears in a typical LJ-like liquid. However, 
the crossover length scale in LJ-like liquids is less sensi- 
tive to temperature and the slope of the dispersed solvo- 
phobes line is greater at lower temperatures, causing the 
region of dcstabilization to dramatically shrink or alto- 
gether disappear. 



Fig. 



10 shows quantitative measures 
of the dissociation size range in the LJ and Jagla liquids 
for cavities equivalent to the solvent size and aggregate 
packing fractions equivalent to the solvent packing frac- 
tion. The dissociation region in the Jagla liquid is orders 
of magnitude larger than that in the LJ liquid. 
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Figure 9. Qualitative depiction of solvation free energy per 
surface area of large solvophobic aggregates and dispersed 
small solutes in (a) typical and (b) water-like solvents. Red 
and blue correspond to warm (Th) and cold temperatures 
(Tl), respectively. The shaded region highlights the aggre- 
gate size range where cooling from Th to Tl destabilizes the 
aggregate. The sloped line which here depicts the rise from 
very small solute to large radius behavior is used to empha- 
size that the shape of this molecular scale transition region is 
represented only generically in this figure. 

In general, the range of the destabilization region is ex- 
tended by cooling to lower temperatures or by composing 
aggregates of smaller constituent particles. A prediction 
made by this model is the possibility of cold-induced dis- 
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Figure 10. The specific case of Fig. [9] for the temperature 
dependence of solvophobic solvation free energies in (a) the 
LJ liquid for T = 0.65 (blue) and T = 0.95 [e L j/k B ] (red) and 
(b) the Jagla liquid for T = 0.4 (blue) and T = 1.0 [e 2 /k B ] 
(red). The constituent solvophobes are equivalent in size to 
the solvent diameter and the aggregate packing fraction is 
taken equivalent to the solvent packing fraction. Both liquids 
have a range of cavity sizes (shaded region) where cooling from 
the warm temperature (red lines) to the cool temperature 
(blue lines) destabilizes the aggregate (solid lines) relative to 
the dispersed spheres (dashed lines). The size range in the 
Jagla liquid is far more pronounced, however (note the order 
of magnitude difference in the abscissa scales). 



sociation of solvophobic aggregates in LJ-like solvents. 
Aggregates composed of sufficiently small cavity solutes 
will in fact, in this model, have a range of sizes for which 
cooling will destabilize the aggregate and induce its de- 
composition. It would indeed be striking if such a limit 
were faithfully captured by this thought experiment in 
spite of its overall simplicity. 



V. CONCLUSIONS 

The results of exhaustive MC simulations of cavity 
formation along the saturation curves of the LJ liquid 
and the Jagla liquid were presented. The temperature- 
dependence of the solvation thermodynamics of cavities 
ranging from one-half to six times the solvent particle 
size were compared between the two simple liquids and 
to predictions for cavity formation in water given by a 
cavity equation of state (C-EoS). The comparisons be- 
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tween the Jagla liquid, water, and the simple liquid (LJ) 
serve to illuminate the features of hydrophobic hydration 
that are unique to water. 

The Jagla liquid demonstrates water-like behavior in 
its resistance to dewetting of large cavity surfaces. In 
the presence of the largest cavity sizes considered (six sol- 
vent diameters), the LJ liquid showed a dewetting transi- 
tion at all thermodynamic states on the saturation curve, 
whereas the Jagla liquid resists dewetting at low temper- 
ature saturation states. 

The Jagla liquid is also water-like in its enthalpic and 
entropic behavior in the sense that the solvation entropy 
of small cavities is negative and the heat capacity incre- 
ment is positive. The LJ liquid on the other hand mani- 
fests a strictly positive entropy for all cavities larger than 
half the solvent size and shows a negligible heat capacity 
increment. 

From our analysis, we infer the important result that 
it is the existence of a second energy scale in the Jagla 
liquid and in water, compared to a simple liquid, that 
energetically favors the creation of void space at low tem- 
peratures, that gives rise to the anomalous liquid state 
properties as well as solvation behavior. Of course, the 
ability of the fluid to access the low energy structures 
with only modest expansion implies that the particular 
length scales involved are closely coupled to this obser- 
vation [9]. 

We have demonstrated that the scaling and tempera- 
ture dependence of the solvation free energies of cavity 
solutes in Jagla liquid is qualitatively similar to that of 
water. Both liquids have negative solvation entropies for 
small cavities that cross over to positive with increasing 



cavity size. These crossovers for the Jagla liquid occur 
at a shorter length scale relative to the solvent size than 
those of water. 

Combining ideas from Chandler [1 and Rajamani et 
al. [5], a simple model for aggregate dissociation was in- 
troduced by modeling an aggregate as a single large hard 
sphere with a volume equal to the sum of the volumes 
of the constituent spheres divided by a packing fraction. 
The consequences of the differing size scaling and tem- 
perature dependence of solvation free energy for the ag- 
gregate compared to the dispersed constituent spheres is 
clearly demonstrated in the context of this simple model 
for aggregation. In particular, it was shown that cold- 
induced dissociation will occur for aggregates composed 
of sufficiently small spheres in water-like liquids. The de- 
gree to which such behavior is accurately described by 
the simple model is of interest for further investigations, 
as is the detailed examination of other two-scale liquids 
containing both a hard and soft core component. 
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